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1.  Introduction.  In  physics  and  engineering  one  often  encounters  systems  de¬ 
scribed  by  differential  equations  which  possess  a  non-trivial  structure  that  embodies 
important  physical  properties  which  can  affect  the  qualitative  behaviour  of  the  sys¬ 
tem.  We  use  the  term  “structure”  in  a  broad  sense  —  we  consider  any  constraint 
upon  the  solutions  as  structure.  Systems  of  ODEs  can  possess  a  wide  range  of  struc¬ 
ture;  for  example,  in  classical  mechanics  the  Poincare  differential  invariants  result  in 
conservation  of  (projections  of)  phase  space  volume.  In  addition,  systems  may  pos¬ 
sess  constants  of  motion  (first  integrals)  such  as  energy  or  angular  momentum,  etc . 
Systems  described  by  PDEs  have  the  extra  complication  that  they  admit  the  possi¬ 
bility  of  an  (uncountable)  infinity  of  invariants  such  as  the  Casimir  invariants  found 
in  continuum  mechanics  [11]-  For  example,  any  functional  of  vorticity  is  conserved  by 
Euler’s  equations,  likewise  any  functional  of  the  phase-space  density  is  conserved  by 
the  Vlasov  equation. 

For  systems  with  structure,  all  numerical  errors  are  not  equal.  One  can  think  of 
the  system’s  structure  as  restricting  the  dynamical  variables  to  some  n-dimensional 
surface  [7]  and  then  separate  numerical  errors  into  two  categories:  those  tangent  to  this 
surface  and  those  normal  to  this  surface.  As  these  errors  accumulate  over  many  time 
steps,  those  in  the  latter  category  are  errors  that  violate  the  structure  of  the  system 
(for  example,  they  may  represent  energy  gain  or  dissipation  in  a  conservative  system) 
while  those  in  the  former  category  are  more  benign  as  they  represent  quantitative 
rather  than  qualitative  errors. 

There  is  much  anecdotal  evidence  (see  for  example  the  discussion  in  Refs.  (8,  13]) 
that  numerical  methods  which  preserve  the  structure  of  a  system  are  likely  to  yield 
results  superior  to  (more  accurate  than)  a  generic  method,  or,  alternatively,  that  for  a 
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given  accuracy,  a  structure  preserving  method  is  likely  to  require  less  computational 
effort.  This  is  in  part  due  to  structure-preserving  methods  having  better  stability 
properties  than  generic  methods;  in  some  cases  structure  preserving  methods  can 
exhibit  unconditional  non-linear  stability.  (Furthermore,  for  the  Korteweg-de  Vries 
equation  it  has  been  rigorously  shown  [8]  that  conservative  methods  exhibit  less  rapid 
error  growth  than  do  non-conserving  methods.)  This  seems  to  be  especially  true  in 
cases  where  the  time  domain  of  interest  is  much  larger  than  the  system’s  characteristic 
time  scale(s). 

These  concepts  are  also  applicable  to  weakly  non-ideal  systems,  i.e.,  those  systems 
that  can  be  viewed  as  a  perturbation  to  an  ideal  system.  For  these  systems  the 
evolution  is  obtained  using  operator  splitting;  one  separates  the  differential  operator 
into  a  piece  describing  the  ideal  system  and  a  piece  describing  the  perturbation.  One 
uses  a  structure-preserving  method  for  the  ideal  part  and  a  generic  method  for  the 
perturbation.  The  complete  evolution  is  obtained  by  combining  the  evolution  for  the 
separate  operators  and  has  the  desirable  property  that  in  the  limit  as  the  perturbation 
vanishes,  the  structure-preserving  behaviour  is  recovered  [4,  14]. 

We  will  proceed  by  considering  several  different  examples  of  structure  preserving 
methods  drawn  from  various  branches  of  physics. 

2.  Symplectic  Integrators.  Consider  a  one-dimensional  Hamiltonian  system: 


(1) 


,  .  dH 

and  p  =  —  . 

dq 


(Throughout,  an  over-dot  will  be  used  to  denote  a  time  derivative.)  Since  this  is  a 
Hamiltonian  system,  phase  space  volume  is  conserved  and  thus  the  Poisson  bracket 
of  q  and  p  at  any  time  is  unity: 


(2) 


[<7(*).P(0] 


_  dg{t)  dp(t) 


dpjt)  dq(t) 
dq (?)  dp(t') 


That  is,  the  dynamical  variables  at  any  one  time  are  related  to  the  dynamical  variables 
at  another  time  by  a  canonical  transformation.  Now,  suppose  one  uses  Euler’s  method 
to  numerically  evolve  this  system.  Denoting  the  numerical  solution  at  time  tn  —  nr 
by  qn  and  pn,  where  the  derivatives  of  H  are  evaluated  at  qn  and  pn, 


(3) 


f)  H 

q"+1=qn-T^(qn,Pn), 

r)H 

pn+l=pn  +  T—(qn,pn)- 

dq 


The  pertinent  question  is  “Does  this  time-advance  map  represent  a  canonical  trans¬ 
formation?”  A  straightforward  calculation  reveals 


(4) 


[qn+\pn+V,P"}  =  1-r2 


d2H 


i  2 


dq  dp\ 


82Hd2H\ 

dq 2  dp 2  J  ’ 


where  again  the  derivatives  of  H  are  evaluated  at  qn  and  pn.  We  see  that  evolution 
given  by  (3)  is  not  Hamiltonian,  Le.,  it  will  not  preserve  phase  space  volume. 

The  solution,  as  was  originally  recognized  by  De  Vogelaere  [9]  in  the  mid  1950’s, 
is  to  make  the  numerical  approximation  to  the  evolution  a  canonical  transformation. 
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The  net  result  is  that  the  numerical  evolution  is  the  exact  solution  of  some  Hamilto¬ 
nian  system  that  approximates  the  original  system.  In  this  way,  the  numerical  error 
normal  to  the  constraint  surface  representing  the  Poincare  invariants  is  eliminated. 
The  importance  of  this  method  for  numerically  evolving  Hamiltonian  systems  was  re¬ 
discovered  in  the  1980’s  by  accelerator  physicists  when  designing  the  Superconducting 
Super  Collider,  where  it  was  necessary  to  understand  the  smearing  of  phase  space1 
due  to  non-linearities  in  the  machine  over  109  particle  orbits.  In  this  case,  symplec- 
tic  methods  were  necessary  both  to  keep  the  computational  demands  reasonable  and 
also  to  ensure  that  the  phase  space  distortions  seen  were  solely  a  consequence  of  the 
machine  design  and  not  numerical  artifacts.  Symplectic  methods  also  figured  promi¬ 
nently  in  the  pioneering  work  of  Wisdom  and  coworkers  in  exploring  the  long-term 
stability  of  the  solar  system,  see  for  example  Ref.  [16]. 

For  a  comprehensive  review  of  these  techniques,  see  the  article  by  Channel  and 
Scovel  [5]  as  well  as  the  monograph  by  Sans-Serna  and  Calvo  [12].  A  large  variety  of 
methods,  both  explicit  (which  are  based  largely  on  operator  splitting  [IT])  and  implicit 
are  now  known.  As  a  simple  example,  the  Mid-Point  rule  (where  the  equations  of 
motion  are  differenced  between  time  steps), 

dH  f  qn+1  +qn  pn+1  +  pn\ 
dp  V  2  ’  2  )  ’ 

dH  f  qn+1  +  qn  pn+1  +  pn\ 
dq  V  2  ’  2  )' 

yields  a  time-advance  map  that  is  a  canonical  transformation. 

3.  Exactly  Conservative  Integrators.  We  now  move  on  to  consider  spectral 
truncations  of  the  Euler  equations.  In  particular,  we  consider  the  “three-wave”  prob¬ 
lem  obtained  by  restricting  the  Fourier-transformed  equations  to  three  modes  [1,  2, 
3,  6,  13]: 

4>k  =  Mk  4>p  4>q  =  Sk{4>)  , 

(6)  <j>p  -  Mp4>q  4>k  =  Sp{4>) , 

4>q  —  Mq  <j>K  4>p  =  Sq(<}>)  , 

where  <j>  =  (<£/c,  </>p,  </>q),  AT,  P,  and  Q  are  the  magnitudes  of  the  Fourier  wavenumbers 
of  the  three  modes,  and  the  coefficients  Mj c,  Afp,  and  Mq  satisfy 


(7) 

Mk  +  Mp  +  Mq  =  0 

and 

(8) 

K2  Mk  +  P2Mp  +  Q2Mq  =  0 

This  system  possesses  two  invariants:  the  total  energy 
(9)  E  =  -  (<j)2K  +  0p  +  <t>o) 

Minimizing  the  phase  space  occupied  by  the  beam  ( “emittance”  in  the  language  of  accelerator 
physicists)  is  critical  to  obtaining  high  luminosity  and  consequently  to  the  usefulness  of  the  collider 
for  physics  experiments. 
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and  the  total  enstrophy 

(10)  Z=l-(K2^K  +  P2<j>2PpQ24>2Q). 

These  invariants  follow  directly  from  properties  of  Sk : 

(Ha)  X)  <t>kSk=0, 

ke{K,p,Q} 

(lib)  P<t>kSk=0. 

ke{K,p,Q} 

(The  equations  (6),  are  identical  to  Euler’s  equations  for  the  rigid  body,  in  which  case 
the  additional  invariant  is  the  norm  of  the  total  angular  momentum.) 

When  (6)  is  integrated  numerically  using  standard  explicit  methods,  neither  E 
nor  Z  are  exactly  conserved.  This  can  be  illustrated  by  applying  Euler’s  method. 
Denoting  the  numerical  solution  at  time  tn  =  nr  by  </?n, 

(12)  V*+1  =  rf  +  rW);  ke{K,P,Q}. 

The  energy  at  the  new  time  is 

E(ipn+1)  =  i  +TSk{<pn)}2 

k 

(13)  =E(^)  +  ^t2X>(^)2, 

k 

where  we  have  used  (11a)  in  the  last  step.  Thus  the  total  energy  is  always  increasing. 
Similarly  we  find 

(14)  Z(^+1)  =  Z(^)+ir2^fc25fc(^)2, 

k 

which  is  likewise  always  increasing. 

Inspired  by  the  idea  of  backwards  error  analysis,  one  is  left  to  wonder  if  it  possible 
to  modify  the  equations  of  motion  such  that,  for  a  given  integrator,  the  time-advance 
mapping  conserves  energy  and  enstrophy  while  still  yielding  a  numerical  approxima¬ 
tion  to  the  original  system  [13].  More  formally,  consider  the  modified  system 

(15)  *b  =  Sife(0)  +  /fc. 

where  fk  is  chosen  to  guarantee  exact  conservation  of  energy  and  enstrophy  and  to 
vanish  in  the  limit  r  — ►  0  sufficiently  rapidly  that  the  numerical  solution  is  still 
consistent  with  the  original  system  (i.e.,  the  order  of  the  new  time-advance  map 
should  be  the  same  as  the  original  integrator). 

We  illustrate  this  analysis  using  a  second-order  predictor-corrector  algorithm 
(Heun’s  method  [18]).  Applying  this  integrator  to  the  modified  system  we  have 

Vk  =  <Pk  +r(Sk(<pn)  +  ffc), 

(16) 

<^fc+1  =  Vk  +  j  (^(V1")  +  f*.4  +  ^fc  +  ffc)  * 
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FlG.  3.1.  Integration  of  the  three-wave  problem  using  a  conventional  second-order  predictor- 
corrector  (dotted  line)  and  the  conservative  predictor-corrector  (solid  line).  Both  methods  took 
approximately  4000  time  steps  of  size  0.05.  Initially  =  \/l75,  pp  —  0,  and  pQ  =  \/1.5.  The 
effect  of  the  4%  energy  gain  by  the  conventional  method  is  clearly  visible.  (From  [13].) 

where  is  the  discretization  of  fk,  Sk  =  Sk((p)  and  =  f k(j?)>  We  expect  that  it 
should  be  sufficient  to  modify  the  corrector: 

(17)  <Pnk+1=‘Pnk  +  l{s*('Pn)  +  Sk+Qk). 

The  conservation  laws  imply 

2  ~  \^Pk  +  2  ( ' Sk(vn )  +  Sk^j\ 

(18)  _ 

+  sgn (<pk)  \J(wl)2  +  t  (pi  Sk(<pn)  +  <pk  Sk )  . 

It  is  straightforward,  if  somewhat  tedious,  to  show  that  gk  —  ®(t2),  thus  preserving 
the  second-order  character  of  the  method.  The  exactly  conservative  time-stepper 
takes  the  form 

<Pk  =  <Pk  +rSk{vn), 

(19) 

VT1  =  sgn(^) 
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Comparing  these  formulas  with  the  original  predictor-corrector  we  see  that  original 
method  is  altered  at  all  orders  in  r  beyond  second. 

In  Fig.  3.1,  we  compare  the  solution  of  (6)  using  the  conventional  predictor-cor¬ 
rector  and  the  conservative  predictor-corrector  for  a  large  number  of  orbits.  Clearly 
the  conservative  method  gives  superior  results.  (See  [13]  for  additional  examples  of 
conservative  integrators  and  for  a  more  thorough  discussion.) 

4.  Unitary  Integrators.  As  a  final  example,  we  consider  an  n-level  quantum 
system  driven  by  external  fields.  The  dynamics  of  the  density  matrix,  p,  generated 
by  a  Hamiltonian  H ,  is  governed  by  the  quantum  Liouville  equation: 

(20)  ihp  =  [H  ,  p] . 

This  equation  has  a  non-trivial  kinematic  structure  —  the  Hioe-Eberly  invariants, 
tr  ffl ,  j  =  1,2,  . . . ,  n,  are  non-evolving  regardless  of  the  form  of  the  Hamiltonian  [10]. 
These  constants  are  a  direct  consequence  of  the  unitary  evolution  of  the  density  matrix 
and  are  the  analogues  of  the  Poincare  invariants  in  classical  mechanics.  A  numerical 
solution  where  these  invariants  are  not  preserved  is  in  danger  of  being  unphysical. 
The  exact  dynamics  proceeds  by  unitary  evolution 

(21)  p(t  +  T)  =  U(t,t  +  T)p{t)U'(t,t  +  T) . 

In  the  spirit  of  symplectic  integrators  discussed  in  §2,  the  kinematic  structure  of  (20) 
will  be  preserved  if  the  numerical  evolution  operator  is  also  unitary. 

We  construct  a  numerical  time-advance  mapping  by  approximating  the  exact 
evolution  while  retaining  the  unitary  property  [15] 

(22)  U(t,  t  +  t)  =  e~ir  Mt'T)  =  U(t,  t  +  r)  +  <D(tk)  , 

where  A  is  a  Hermitian  matrix  that  will  depend  on  the  Hamiltonian. 

By  matching  the  Taylor  series  solution  of  (20)  term-by-term  in  r  with  the  ap¬ 
proximate  evolution  generated  by  [/(£,£  4-  r),  the  following  approximations  for  A(t,r) 
are  obtained:  to  second  order 

(23)  A  =  H  +  t  i/; 


to  third  order 


(24)  +  +  + 

and  to  fourth  order 


(25)  -4-ff  +  5rff  +  5r ■f'  +  jr 


where  [•  ,  •]  is  the  matrix  commutator.  Note  that  to  obtain  accuracy  beyond  second- 
order,  one  must  take  into  account  that,  in  general,  [H(ti)  ,  H(t2)}  ^  0. 

Fortunately,  to  use  these  expressions  for  U,  it  is  not  necessary  to  exponentiate  a 
(general)  nxn  matrix.  We  are  free  to  approximate  A  in  any  way  consistent  with  the 
order  of  the  method.  As  described  in  [14]  and  [15]  it  is  possible  to  introduce  a  suitable 
basis,  {Ajt}^!,  such  that  exponentials  of  the  basis  matrices  are  easily  computed.  We 
then  write 


(26) 


e-irA(t,T)  __ 


k~\ 
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where  the  7*  are  determined  by  matching  terms  on  each  side  of  (26)  order-by-order 
in  r  through  G(rK~2). 

As  an  example,  consider  a  simple  two-level  system  with  a  time-independent  Hamil¬ 
tonian 


(27) 


A  second-order  integrator  is  given  by 


(28) 
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In  Fig.  4.1(a)  and  (6),  we  show  the  numerical  solution  of  (20)  obtained  using  the 
unitary  integrator  (28)  as  well  as  that  obtained  with  a  second-order  predictor-correc¬ 
tor.  The  parameters  for  this  example  are  e  =  1,  u  =  0.01,  and  r  =  0.1  and  the  initial 
condition  is 

1(1  e-i7r/4 

(29)  P(0)=2  (eW4  %  J 


In  Fig.  4.1(c)  and  (d)  we  show  the  difference  between  the  solutions  shown  in  (a)  and  (6) 
and  a  high-accuracy  integration.  It  is  evident  that  the  error  of  the  unitary  integrator 
is  much  smaller  than  that  of  the  predictor-corrector  even  though  both  methods  used 
the  same  time  step  and  are  both  second  order.  While  both  methods  conserve  tr  p  (it 
is  a  linear  invariant,  and  hence  preserved  by  predictor-corrector),  this  alone  is  not 
sufficient  to  reproduce  the  dynamics  faithfully. 
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FlG.  4.1.  Results  of  solving  (20)  using  the  unitary  integrator  (solid  line)  and  a  second-order 
predictor-corrector  (heavy  dashes):  numerical  solution  pn  (a)  and  Im  pX2  (6);  numerical  errors 
in  pi i  (c)  and  Im  p\2  ( d ). 


5.  Conclusion.  For  systems  of  differential  equations  that  possess  structure,  bet¬ 
ter  numerical  results  are  obtained  when  the  numerical  algorithms  used  to  solve  these 
equations  respect  this  structure.  For  systems  with  non-linear  algebraic  invariants,  we 
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have  shown  that  it  is  possible  to  systematically  construct  explicit  exactly  conservative 
algorithms.  Furthermore,  we  have  seen  that  kinematic  structure  is  preserved  when  the 
numerical  time-advance  map  shares  the  group  properties  of  the  exact  dynamics:  for 
Hamiltonian  mechanics,  we  preserve  the  Poincare  invariants  when  the  time-advance 
is  a  canonical  transformation;  for  the  quantum  Liouville  equation,  the  Hioe-Eberly 
invariants  are  preserved  when  the  time-advance  is  a  unitary  transformation.  Not  only 
are  the  numerical  results  superior,  but  by  using  structure  preserving  methods,  we  are 
guaranteed  the  remaining  numerical  errors  do  not  violate  the  inherent  physics  of  the 
models.  The  key  point  is  that  integrators  should  have  maximum  knowledge  of  the 
system  they  are  being  used  to  solve . 
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LABORATORY  OPERATIONS 


The  Aerospace  Corporation  functions  as  an  “architect-engineer”  for  national  security  programs,  specializing  in 
advanced  military  space  systems.  The  Corporation's  Laboratory  Operations  supports  the  effective  and  timely 
development  and  operation  of  national  security  systems  through  scientific  research  and  the  application  of 
advanced  technology.  Vital  to  the  success  of  the  Corporation  is  the  technical  staffs  wide-ranging  expertise  and 
its  ability  to  stay  abreast  of  new  technological  developments  and  program  support  issues  associated  with  rapidly 
evolving  space  systems.  Contributing  capabilities  are  provided  by  these  individual  organizations: 

Electronics  and  Photonics  Laboratory:  Microelectronics,  VLSI  reliability,  failure 
analysis,  solid-state  device  physics,  compound  semiconductors,  radiation  effects,  infrared 
and  CCD  detector  devices,  data  storage  and  display  technologies;  lasers  and  electro-optics, 
solid  state  laser  design,  micro-optics,  optical  communications,  and  fiber  optic  sensors; 
atomic  frequency  standards,  applied  laser  spectroscopy,  laser  chemistry,  atmospheric 
propagation  and  beam  control,  LIDAR/LADAR  remote  sensing;  solar  cell  and  array  testing 
and  evaluation,  battery  electrochemistry,  battery  testing  and  evaluation. 

Space  Materials  Laboratory:  Evaluation  and  characterizations  of  new  materials  and 
processing  techniques:  metals,  alloys,  ceramics,  polymers,  thin  films,  and  composites; 
development  of  advanced  deposition  processes;  nondestructive  evaluation,  component 
failure  analysis  and  reliability;  structural  mechanics,  fracture  mechanics,  and  stress 
corrosion;  analysis  and  evaluation  of  materials  at  cryogenic  and  elevated  temperatures; 
launch  vehicle  fluid  mechanics,  heat  transfer  and  flight  dynamics;  aerothermodynamics; 
chemical  and  electric  propulsion;  environmental  chemistry;  combustion  processes;  space 
environment  effects  on  materials,  hardening  and  vulnerability  assessment;  contamination, 
thermal  and  structural  control;  lubrication  and  surface  phenomena. 

Space  Science  Application  Laboratory:  Magnetospheric,  auroral  and  cosmic  ray  physics, 
wave-particle  interactions,  magnetospheric  plasma  waves;  atmospheric  and  ionospheric 
physics,  density  and  composition  of  the  upper  atmosphere,  remote  sensing  using 
atmospheric  radiation;  solar  physics,  infrared  astronomy,  infrared  signature  analysis; 
infrared  surveillance,  imaging,  remote  sensing,  and  hyperspectral  imaging;  effects  of  solar 
activity,  magnetic  storms  and  nuclear  explosions  on  the  Earth’s  atmosphere,  ionosphere  and 
magnetosphere;  effects  of  electromagnetic  and  particulate  radiations  on  space  systems; 
space  instrumentation,  design  fabrication  and  test;  environmental  chemistry,  trace  detection; 
atmospheric  chemical  reactions,  atmospheric  optics,  light  scattering,  state-specific  chemical 
reactions  and  radiative  signatures  of  missile  plumes. 

Center  for  Microtechnology:  Microelectromechanical  systems  (MEMS)  for  space 
applications;  assessment  of  microtechnology  space  applications;  laser  micromachining; 
laser-surface  physical  and  chemical  interactions;  micropropulsion;  micro-  and 
nanosatellite  mission  analysis;  intelligent  microinstruments  for  monitoring  space  and 
launch  system  environments. 

Office  of  Spectral  Applications:  Multispectral  and  hyperspectral  sensor  development; 
data  analysis  and  algorithm  development;  applications  of  multispectral  and  hyperspectral 
imagery  to  defense,  civil  space,  commercial,  and  environmental  missions. 


